
## Playing Politics with Environmental Protection: The Political Economy of Designating Protected Areas
## Jorge Mangonnet, Jacob Kopas, and Johannes Urpelainen
## August 2021

## File: code for producing regression tables (.tex format) and plots 
## NOTE: first run "analysis.R" script to prepare data


### Preamble -------------------------------------------------------------------

## Load packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  maptools,
  spatstat,
  rgdal,
  sp,
  spdep,
  rgeos,
  foreign,
  stargazer,
  xtable,
  ggplot2,
  plyr,
  dplyr,
  reshape2,
  lfe,
  plm,
  stringr,
  scales,
  ggExtra,
  magrittr,
  coefplot
  )

## Functions for regression models and plots
source(paste0(user, "manuscript/replication/Brazil_func.R"))



### (A) TABLES: MAIN MODELS ----------------------------------------------------

## Table 1: Federal Protected Areas and Political Alignment ----
dep.labels <- "Federal Protected Area"
cov.labs <- c("Coalition Alignment", "Party Alignment", "Fed. Prot. Area ('97)")
col.sep <- c(3,3)
col.labels <- c("Federal Protected Area", "Federal Protected Area")
sink(paste0("tables/coal.tex"))
out.tab(c(coal, party), modeltype = 0, c(coal.ug, party.ug))
sink()


## Table 2: Presidential Vote Share and Federal Protected Areas ----
dep.labels <- NULL
cov.labs <- NULL
col.sep <- NULL
cov.labs <- c("Fed. Prot. Area", "Fed. Prot. Area ('97)", 
              "Coalition Alignment", "Fed. Prot. Area:Alignment")
omit.list <- "Constant"
col.sep <- c(3,3)
col.labels <- c("President Party Vote Share", "President Party Vote Share")
sink(paste0("tables/pres.tex"))
out.tab(pres, modeltype = 8, pres.ug)
sink()



### (B) TABLES: EXPLORATORY ANALYSES -------------------------------------------

## Table A5: Comparison of Extensive and Intensive Margins of Coalition Alignment ----
dep.labels <- NULL
cov.labs <- c("Coalition Alignment", "Dummy Fed. Prot. Area ('97)", "Prop. Fed. Prot. Area ('97)")
col.sep <- c(3, 3)
col.labels <- c("Federal Protected Area (Dummy)", "Federal Protected Area (Proportion)")
sink(paste0("tables/muni.tex"))
out.tab(muni, modeltype = 5, muni.ug)
sink()


## Table A6: Interaction between Coalition Alignment and Presidential Vote Share ----
dep.labels <- NULL
cov.labs <- c("Coalition Alignment", "Fed. Prot. Area ('97)", 
              "Pres. Vote Share", "Alignment:Vote Share")
col.sep <- 3
col.labels <- c("Federal Protected Area")
sink(paste0("tables/core_swing.tex"))
out.tab(core, modeltype = 1, core.ug)
sink()


## Table A7: Interaction between Coalition Alignment and Margin of Victory ----
## Targeting of core/swing districts using vote margin
dep.labels <- NULL
cov.labs <- c("Coalition Alignment", "Fed. Prot. Area ('97)", 
              "Pres. Vote Margin", "Alignment:Vote Margin")
col.sep <- 3
col.labels <- c("Federal Protected Area")
sink(paste0("tables/core_swing_margin.tex"))
out.tab(vmarg, modeltype = 1, vmarg.ug)
sink()


## Table A8: Difference-in-Differences Estimation of the Effect of Protected Areas on Local Extractive Industries ----
dep.labels <- NULL
cov.labs <- c("Protected Area ('97)", "Post 2001", "Protected Area ('97):Post 2001")
omit.list <- "Constant"
col.sep <- c(3,3)
col.labels <- c("Soybean Production", "Mining Leases")
sink(paste0("tables/muni_extraction.tex"))
out.tab(ex, modeltype = 9, ex.ug)
sink()


## Table A9: Interaction between Coalition Alignment and Potential for Economic Exploitation ----
dep.labels <- NULL
cov.labs <-c("Coalition Alignment", "Fed. Prot. Area ('97)", 
             "Deforested", "Alignment:Deforested", 
             "Pastures", "Alignment:Pastures", 
             "Soybean", "Alignment:Soybean")
col.sep <- 3
col.labels <- c("Federal Protected Area")
sink(paste0("tables/interaction.tex"))
out.tab(int, modeltype = 3, int.ug)
sink()


## Table A10: Interaction between Coalition Alignment and Federal Roadways ----
dep.labels <- NULL
cov.labs <-c("Coalition Alignment", "Fed. Prot. Area ('97)", 
             "Transamaz\\^{o}nica", "Alignment:Transamaz\\^{o}nica", 
             "Main Federal Highway", "Alignment:Main Federal Highway", 
             "Federal Road", "Alignment:Federal Road")
col.sep <- 3
col.labels <- c("Federal Protected Area", "Federal Protected Area")
sink(paste0("tables/interaction_rds.tex"))
out.tab(int.rds, modeltype = 4, int.rds.ug)
sink()



### (C) TABLES: PLACEBO TESTS AND ROBUSTNESS CHECKS ----------------------------

## Table A11: Effect of Coalition Alignment on Federal Protected Areas, State Protected Areas, and Indigenous Lands ----
dep.labels <- NULL
cov.labs <- c("Coalition Alignment", "Fed. Prot. Area ('97)", "Indigenous Lands ('97)", "State Prot. Area ('97)")
col.sep <- c(3,3,3)
col.labels <- c("Federal Protected Area", "Indigenous Lands", "State Protected Area")
sink(paste0("tables/type.tex"))
out.tab(type, modeltype = 2, type.ug)
sink()


## Table A12: State Protected Areas and Governor-Mayor Alignment ----
dep.labels <- "State Protected Area"
cov.labs <- c("State Coalition Alignment", "State Party Alignment", "State. Prot. Area ('97)")
col.sep <- c(3,3)
col.labels <- c("State Protected Area", "State Protected Area")
sink(paste0("tables/state.tex"))
out.tab(c(s_coal, s_party), modeltype = 0, c(scoal.ug, sparty.ug))
sink()


## Table A13: Controlling for Distance to the Border ----
dep.labels <- NULL
cov.labs <- c("Coalition Alignment", "Fed. Prot. Area ('97)", 
              "Distance", "Distance:Alignment")
col.sep <- c(2)
col.labels <- c("Federal Protected Area")
sink(paste0("tables/rd.tex"))
out.tab(rd, modeltype = 4, cxrd.ug)
sink()


## Table A14: Controlling for Geographic Coordinates ----
dep.labels <- NULL
cov.labs <- c("Coalition Alignment", "Fed. Prot. Area ('97)", 
              "x-coord.", "y-coord.", 
              "x$^2$", "y$^2$", 
              "x $\\cdot$ y", 
              "x$^3$", "y$^3$", 
              "x$^2 \\cdot$ y", 
              "x $\\cdot$ y$^2$")
omit.list <- "Constant"
col.sep <- 4
col.labels <- c("Federal Protected Area")
sink(paste0("tables/dell.tex"))
out.tab(dell, modeltype = 11, dell.ug)
sink()


## Table A15: Nonparametric Estimation of the Effect of Coalition Alignment ----
## NOTE: stargazer function incompatible w/ rdrobust;
##       first create simple lm object and then replace estimates w/ those of rdrobust
# Create lm object for storing estimates
nplm0 <- lm(federal ~ p_coal, data = df.s)
# Get coefficients, standard errors, z-scores, and p-values
# [1] Conventional, [2] Bias-Corrected, [3] Robust
np.coef.c <- np$coef[1]
np.coef.b <- np$coef[2]
np.coef.r <- np$coef[3]
np.se.c <- np$se[1]
np.se.b <- np$se[2]
np.se.r <- np$se[3]
np.z.c <- np$z[1]
np.z.b <- np$z[2]
np.z.r <- np$z[3]
np.pv.c <- np$pv[1]
np.pv.b <- np$pv[2]
np.pv.r <- np$pv[3]
# Table arguments
dep.labels <- NULL
omit.list <- NULL
col.sep <- c(1,1,1)
cov.labs <- c("Coalition Alignment")
col.labels <- c("\\textit{Conventional}", "\\textit{Bias-Corrected}", "\\textit{Robust}")
dep.labels <- "Federal Protected Area"
add.lines <- list(c("Bandwidth", 25, 25, 25),
                  c("Lower C.I.", round(np$ci[1], 3), round(np$ci[2], 3), round(np$ci[3], 3)),
                  c("Upper C.I.", round(np$ci[4], 3), round(np$ci[5], 3), round(np$ci[6], 3)),
                  c("Treated obs.", np$N[1], np$N[1], np$N[1]),
                  c("Control obs.", np$N[2], np$N[2], np$N[2]))
stargazer(nplm0, nplm0, nplm0,
          coef = list(c(0, np.coef.c), c(0, np.coef.b), c(0, np.coef.r)),
          se = list(c(0, np.se.c), c(0, np.se.b), c(0, np.se.r)),
          t = list(c(0, np.z.c), c(0, np.z.b), c(0, np.z.r)),
          p = list(c(0, np.pv.c), c(0, np.pv.b), c(0, np.pv.r)),
          header = F,
          float = F,
          no.space = T,
          covariate.labels = cov.labs,
          keep = "p_coal",
          omit.stat = c("rsq", "adj.rsq", "ser", "f", "n"),
          add.lines = add.lines,
          column.labels = col.labels,
          dep.var.labels = dep.labels,
          dep.var.labels.include = TRUE,
          star.cutoffs = c(0.10, 0.05, 0.01, 0.001),
          star.char = c("+", "*", "**", "***"),
          notes = c("$^{+}$p<0.1; $^{*}$p<0.05; $^{**}$p<0.01; $^{***}$p<0.001"),
          notes.append = FALSE,
          column.separate = col.sep, 
          out = paste0("tables/np.tex"))


## Table A16: Comparison of Main Models at Different Bandwidths ----
dep.labels <- "Federal Protected Area"
cov.labs <- c("Coalition Alignment", "Fed. Prot. Area ('97)")
col.sep <- c(3, 3, 3)
col.labels <- c("20km Bandwidth", "15km Bandwidth", "10km Bandwidth")
sink(paste0("tables/band.tex"))
out.tab(bands, modeltype = 2, bands.ug)
sink()


## Table A17: Main Models with Full Panel of Grid Cell-Year Observations ----
dep.labels <- "Federal Protected Area"
cov.labs <- c("Coalition Alignment", "Fed. Prot. Area ('97)")
col.sep <- 3
col.labels <- "Federal Protected Area"
sink(paste0("tables/full.tex"))
out.tab(full, modeltype = 1, full.ug)
sink()


## Table A18: Main Models without Fully Saturated Grid Cells ----
dep.labels <- "Federal Protected Area"
cov.labs <- c("Coalition Alignment", "Fed. Prot. Area ('97)")
col.sep <- 3
col.labels <- "Federal Protected Area"
sink(paste0("tables/satur.tex"))
out.tab(sat, modeltype = 1, sat.ug)
sink()


## Table A19: Main Models with Imbalanced Pre-treatment Covariates as Control Variables ----
## NOTE: pre-treatment covariates not shown in the final table ('omit.list' argument)
dep.labels <- "Federal Protected Area"
cov.labs <- c("Coalition Alignment", "Fed. Prot. Area ('97)")
omit.list <- c("ideo_imp_94", "deforested", "urb_extent", "sgi_mean", "elf", "cfr_mean", "acc_mean", "rcr_mean", "sbi_mean", "sbr_mean", "mml_mean", "cfi_mean", "rci_mean")
col.sep <- c(2,2)
col.labels <- "Federal Protected Area"
sink(paste0("tables/coal_cv.tex"))
out.tab(ccv, modeltype = 12, ccv.ug)
sink()



### (D) TABLES: EXTENSIONS -------------------------------------------------------

## Table A20: Environmental Embargoes and Political Alignment ----
dep.labels <- "Environmental Embargoes"
cov.labs <- c("Coalition Alignment", "Party Alignment", "Fed. Prot. Area ('97)")
col.sep <- c(3,3)
col.labels <- c("Environmental Embargoes", "Environmental Embargoes")
sink(paste0("tables/embarg.tex"))
out.tab(c(embc, embp), modeltype = 0, c(embc.ug, embp.ug))
sink()


## Table A21: Federal Protected Areas and Coalition Alignment in the PSDB Presidential Term (1997-2002) ----
dep.labels <- NULL
cov.labs <- c("Coalition Alignment", "Fed. Prot. Area ('97)")
col.sep <- 3
col.labels <- "Fed. Protected Area"
sink(paste0("tables/cardoso.tex"))
out.tab(cardoso, modeltype = 1, card.ug)
sink()


## Table A22: Federal Protected Areas and Coalition Alignment in the PT Presidential Terms (2003-2012) ----
dep.labels <- NULL
cov.labs <- c("Coalition Alignment", "Fed. Prot. Area ('97)")
col.sep <- 3
col.labels <- "Fed. Protected Area"
sink(paste0("tables/lula.tex"))
out.tab(lula, modeltype = 1, lula.ug)
sink()


## Table A23: Incumbent Mayor’s Vote Share and Federal Protected Areas ----
dep.labels <- NULL
cov.labs <- c("Fed. Prot. Area", "Fed. Prot. Area ('97)", "Coalition Alignment", "Alignment:Prot. Area.")
omit.list <- "Constant"
col.sep <- 6
col.labels <- "Mayor's Vote Share"
sink(paste0("tables/mayor.tex"))
out.tab(mayor, modeltype = 9, mayor.ug)
sink()


## Table A24: Difference-in-Differences Estimation of the Effect of Protected Areas on Local Peasant Agriculture ----
dep.labels <- NULL
cov.labs <- c("Protected Area ('97)", "Post 2001", "Protected Area ('97):Post 2001")
omit.list <- "Constant"
col.sep <- c(3,3)
col.labels <- c("Peasant Families", "Peasant Farmsteads")
sink(paste0("tables/muni_peasant.tex"))
out.tab(pea, modeltype = 9, pea.ug)
sink()


## Table A25: Interaction between Coalition Alignment and Municipal Deforestation ----
dep.labels <- NULL
cov.labs <- c("Coalition Alignment", "Fed. Prot. Area ('97)", 
              "Deforested Municipal", "Alignment:Deforested Municipal", 
              "Critical Area", "Alignment:Critical Area")
col.sep <- c(3,3)
col.labels <- c("Federal Protected Area", "Federal Protected Area")
sink(paste0("tables/md.tex"))
out.tab(md, modeltype = 0, md.ug)
sink()


## Table A26: Federal Protected Areas and President-Governor Alignment ----
dep.labels <- "Federal Protected Area"
cov.labs <- c("President-Governor", "Fed. Prot. Area ('97)")
col.sep <- 3
col.labels <- "Federal Protected Area"
sink(paste0("tables/presgov.tex"))
out.tab(g_coal, modeltype = 10, gcoal.ug)
sink()


### (E) OTHER TABLES -----------------------------------------------------------
### Sample Construction
### Identifying Assumptions

## Table A1: Summary Statistics ----
sum.tab <- data.frame(years = paste0(min(df.s$year), " - ", max(df.s$year)),
                      obs = nrow(df.s),
                      grids = length(unique(df.s$ID)),
                      munis = length(unique(df.s$codigo)),
                      muni.p = length(unique(df.s$m_pair)),
                      prot = length(which(df.s$federal > 0)),
                      mprot = round(mean(df.s$federal), digits = 4),
                      tobs = nrow(df.s[which(df.s$p_coal == 1), ]),
                      tprot = length(which(df.s$federal[which(df.s$p_coal == 1)] > 0)),
                      tmprot = round(mean(df.s$federal[df.s$p_coal ==1]), digits = 4),
                      cobs = nrow(df.s[which(df.s$p_coal == 0), ]),
                      cprot = length(which(df.s$federal[which(df.s$p_coal == 0)] > 0)),
                      cmprot = round(mean(df.s$federal[df.s$p_coal == 0]), digits = 4))
sum.tab <- t(sum.tab)
sum.tab[2:nrow(sum.tab), ] <- prettyNum(sum.tab[2:nrow(sum.tab), ], big.mark = ",")
colnames(sum.tab) <- "Coalition Alignment"
row.names(sum.tab) <- c("Years", "Observations", "Grid Cells", "Municipalities", "Muni. Pairs", "Protected Area", "Mean Prot. Area",
                        "Observations ","Protected Area ", "Mean Prot. Area ",
                        "Observations  ","Protected Area  ", "Mean Prot. Area  ")

addtorow <- list()
addtorow$pos <- list( 1, 7, 10)
addtorow$command <- c("\\\\\n \\emph{All Observations} \\\\\n",
                      "\\\\\n \\emph{Treated Observations} \\\\\n",
                      "\\\\\n \\emph{Control Observations} \\\\\n")

sink("tables/sum_stat.tex")
print(xtable(sum.tab, auto = TRUE),
      include.rownames = TRUE, include.colnames = FALSE, 
      add.to.row = addtorow,
      floating = FALSE, booktabs = TRUE)
sink()


## Table A3: Balance Test for Coalition Alignment on 40 Pre-Treatment Covariates in the 2000, 2004, and 2008 Mayoral Elections ----
cb.tab <- data.frame(
  # Difference in means (2000)
  diff00=c(c00.ideol$coefficients[1], c00.defo$coefficients[1], c00.rain$coefficients[1], c00.eva$coefficients[1], c00.temp$coefficients[1],
              c00.climaggr$coefficients[1], c00.alti$coefficients[1], c00.slope$coefficients[1], 
              c00.access$coefficients[1], c00.workab$coefficients[1], c00.nutri$coefficients[1], 
              c00.water$coefficients[1], c00.veget$coefficients[1], c00.cacao_i$coefficients[1], 
              c00.cacao_r$coefficients[1], c00.coffee_i$coefficients[1], c00.coffee_r$coefficients[1], 
              c00.grass_i$coefficients[1], c00.grass_r$coefficients[1], c00.legum_i$coefficients[1], 
              c00.legum_r$coefficients[1], c00.maize_i$coefficients[1], c00.maize_r$coefficients[1], 
              c00.rice_i$coefficients[1], c00.rice_r$coefficients[1], c00.soy_i$coefficients[1], 
              c00.soy_r$coefficients[1], c00.sugar_i$coefficients[1], c00.sugar_r$coefficients[1], 
              c00.amphib$coefficients[1], c00.birds$coefficients[1], c00.mammal$coefficients[1],
              c00.urbdist$coefficients[1], c00.urbexte$coefficients[1], 
              c00.br230$coefficients[1], c00.mainrds$coefficients[1], c00.roads$coefficients[1], 
              c00.elf$coefficients[1], c00.popc$coefficients[1], c00.popd$coefficients[1]),
  # P-values (2000)
  pvalue_00=c(c00.ideol$cpval, c00.defo$cpval, c00.rain$cpval, c00.eva$cpval, c00.temp$cpval,
                 c00.climaggr$cpval, c00.alti$cpval, c00.slope$cpval, 
                 c00.access$cpval, c00.workab$cpval, c00.nutri$cpval, 
                 c00.water$cpval, c00.veget$cpval, c00.cacao_i$cpval, 
                 c00.cacao_r$cpval, c00.coffee_i$cpval, c00.coffee_r$cpval, 
                 c00.grass_i$cpval, c00.grass_r$cpval, c00.legum_i$cpval, 
                 c00.legum_r$cpval, c00.maize_i$cpval, c00.maize_r$cpval, 
                 c00.rice_i$cpval, c00.rice_r$cpval, c00.soy_i$cpval, 
                 c00.soy_r$cpval, c00.sugar_i$cpval, c00.sugar_r$cpval, 
                 c00.amphib$cpval, c00.birds$cpval, c00.mammal$cpval, 
                 c00.urbdist$cpval, c00.urbexte$cpval, 
                 c00.br230$cpval, c00.mainrds$cpval, c00.roads$cpval, 
                 c00.elf$cpval, c00.popc$cpval, c00.popd$cpval),
  # Difference in means (2004)
  diff04=c(c04.ideol$coefficients[1], c04.defo$coefficients[1], c04.rain$coefficients[1], c04.eva$coefficients[1], c04.temp$coefficients[1],
              c04.climaggr$coefficients[1], c04.alti$coefficients[1], c04.slope$coefficients[1], 
              c04.access$coefficients[1], c04.workab$coefficients[1], c04.nutri$coefficients[1], 
              c04.water$coefficients[1], c04.veget$coefficients[1], c04.cacao_i$coefficients[1], 
              c04.cacao_r$coefficients[1], c04.coffee_i$coefficients[1], c04.coffee_r$coefficients[1], 
              c04.grass_i$coefficients[1], c04.grass_r$coefficients[1], c04.legum_i$coefficients[1], 
              c04.legum_r$coefficients[1], c04.maize_i$coefficients[1], c04.maize_r$coefficients[1], 
              c04.rice_i$coefficients[1], c04.rice_r$coefficients[1], c04.soy_i$coefficients[1], 
              c04.soy_r$coefficients[1], c04.sugar_i$coefficients[1], c04.sugar_r$coefficients[1], 
              c04.amphib$coefficients[1], c04.birds$coefficients[1], c04.mammal$coefficients[1],
              c04.urbdist$coefficients[1], c04.urbexte$coefficients[1], 
              c04.br230$coefficients[1], c04.mainrds$coefficients[1], c04.roads$coefficients[1], 
              c04.elf$coefficients[1], c04.popc$coefficients[1], c04.popd$coefficients[1]),
  # P-values (2004)
  pvalue04=c(c04.ideol$cpval, c04.defo$cpval, c04.rain$cpval, c04.eva$cpval, c04.temp$cpval,
                 c04.climaggr$cpval, c04.alti$cpval, c04.slope$cpval, 
                 c04.access$cpval, c04.workab$cpval, c04.nutri$cpval, 
                 c04.water$cpval, c04.veget$cpval, c04.cacao_i$cpval, 
                 c04.cacao_r$cpval, c04.coffee_i$cpval, c04.coffee_r$cpval, 
                 c04.grass_i$cpval, c04.grass_r$cpval, c04.legum_i$cpval, 
                 c04.legum_r$cpval, c04.maize_i$cpval, c04.maize_r$cpval, 
                 c04.rice_i$cpval, c04.rice_r$cpval, c04.soy_i$cpval, 
                 c04.soy_r$cpval, c04.sugar_i$cpval, c04.sugar_r$cpval, 
                 c04.amphib$cpval, c04.birds$cpval, c04.mammal$cpval, 
                 c04.urbdist$cpval, c04.urbexte$cpval, 
                 c04.br230$cpval, c04.mainrds$cpval, c04.roads$cpval, 
                 c04.elf$cpval, c04.popc$cpval, c04.popd$cpval),
  # Difference in means (2008)
  diff08=c(c08.ideol$coefficients[1], c08.defo$coefficients[1], c08.rain$coefficients[1], c08.eva$coefficients[1], c08.temp$coefficients[1],
              c08.climaggr$coefficients[1], c08.alti$coefficients[1], c08.slope$coefficients[1], 
              c08.access$coefficients[1], c08.workab$coefficients[1], c08.nutri$coefficients[1], 
              c08.water$coefficients[1], c08.veget$coefficients[1], c08.cacao_i$coefficients[1], 
              c08.cacao_r$coefficients[1], c08.coffee_i$coefficients[1], c08.coffee_r$coefficients[1], 
              c08.grass_i$coefficients[1], c08.grass_r$coefficients[1], c08.legum_i$coefficients[1], 
              c08.legum_r$coefficients[1], c08.maize_i$coefficients[1], c08.maize_r$coefficients[1], 
              c08.rice_i$coefficients[1], c08.rice_r$coefficients[1], c08.soy_i$coefficients[1], 
              c08.soy_r$coefficients[1], c08.sugar_i$coefficients[1], c08.sugar_r$coefficients[1], 
              c08.amphib$coefficients[1], c08.birds$coefficients[1], c08.mammal$coefficients[1],
              c08.urbdist$coefficients[1], c08.urbexte$coefficients[1], 
              c08.br230$coefficients[1], c08.mainrds$coefficients[1], c08.roads$coefficients[1], 
              c08.elf$coefficients[1], c08.popc$coefficients[1], c08.popd$coefficients[1]),
  # P-values (2008)
  pvalue08=c(c08.ideol$cpval, c08.defo$cpval, c08.rain$cpval, c08.eva$cpval, c08.temp$cpval,
                 c08.climaggr$cpval, c08.alti$cpval, c08.slope$cpval, 
                 c08.access$cpval, c08.workab$cpval, c08.nutri$cpval, 
                 c08.water$cpval, c08.veget$cpval, c08.cacao_i$cpval, 
                 c08.cacao_r$cpval, c08.coffee_i$cpval, c08.coffee_r$cpval, 
                 c08.grass_i$cpval, c08.grass_r$cpval, c08.legum_i$cpval, 
                 c08.legum_r$cpval, c08.maize_i$cpval, c08.maize_r$cpval, 
                 c08.rice_i$cpval, c08.rice_r$cpval, c08.soy_i$cpval, 
                 c08.soy_r$cpval, c08.sugar_i$cpval, c08.sugar_r$cpval, 
                 c08.amphib$cpval, c08.birds$cpval, c08.mammal$cpval, 
                 c08.urbdist$cpval, c08.urbexte$cpval,
                 c08.br230$cpval, c08.mainrds$cpval, c08.roads$cpval,
                 c08.elf$cpval, c08.popc$cpval, c08.popd$cpval), 
  # Covariate labels
  row.names=c("Municipal ideological score", "Deforested area", "Rainfall", "Evapotranspiration", "Temperature",
              "Climate aggressiveness", "Altitude", "Slope",
              "Accessibility", "Workability", "Nutrients",
              "Water bodies", "Vegetation", "Cacao suit. (irrigated)",
              "Cacao suit. (rain-fed)", "Coffee suit. (irrigated)", "Coffee suit. (rain-fed)", 
              "Pasture grasses suit. (irrigated)", "Pasture grasses suit. (rain-fed)", "Pasture legumes suit. (irrigated)",
              "Pasture legumes suit. (rain-fed)", "Maize suit. (irrigated)", "Maize suit. (rain-fed)",
              "Rice suit. (irrigated)",  "Rice suit. (rain-fed)", "Soybeans suit. (irrigated)",
              "Soybeans suit. (rain-fed)", "Sugarcane suit. (irrigated)", "Sugarcane suit. (rain-fed)",
              "Threatened amphibians", "Threatened birds", "Threatened mammals",
              "Min. dist. from urban areas", "Urban", 
              "Min. dist. from Transamazonica", "Min. dist. from federal highway", "Min. dist. from federal road",
              "Ethnolinguistic fractionalization index", "Population count", "Population density")
)

addtorow <- list()
addtorow$pos <- list(0, 0, 40)
addtorow$command <- c(paste0(paste0("& \\multicolumn{2}{c}{2000} & \\multicolumn{2}{c}{2004} & \\multicolumn{2}{c}{2008}", collapse=''), '\\\\'),
                      paste0(paste0("& diff. & p-val. & diff. & p-val. & diff. & p-val.", collapse=''), '\\\\'),
                      paste0(paste0(" \\emph{N} & \\multicolumn{2}{c}{52918} & \\multicolumn{2}{c}{64194} & \\multicolumn{2}{c}{64664}", collapse=''), '\\\\')
                      )

sink("tables/bal_stats.tex")
print(xtable(cb.tab, align=c("l", rep("c",6)), digits=3),
      include.colnames = FALSE, 
      add.to.row = addtorow, 
      floating = FALSE, 
      booktabs = TRUE,
      size = "small",
      table.placement="h!"
)
sink()


## Table A4: Moran’s I Test for Coalition Alignment on 40 Pre-Treatment Covariates in the 1996, 2000, 2004, and 2008 Mayoral Elections ----
mi.tab <- data.frame(
  # Moran's I (1996)
  test96=c(round(mi96.ideol$estimate[[1]], 3),
           round(mi96.defo$estimate[[1]], 3), round(mi96.rain$estimate[[1]], 3), 
           round(mi96.eva$estimate[[1]], 3), round(mi96.temp$estimate[[1]], 3),
           round(mi96.climaggr$estimate[[1]], 3), round(mi96.alti$estimate[[1]], 3), round(mi96.slope$estimate[[1]], 3), 
           round(mi96.access$estimate[[1]], 3), round(mi96.workab$estimate[[1]], 3), round(mi96.nutri$estimate[[1]], 3), 
           round(mi96.water$estimate[[1]], 3), round(mi96.veget$estimate[[1]], 3), round(mi96.cacao_i$estimate[[1]], 3), 
           round(mi96.cacao_r$estimate[[1]], 3), round(mi96.coffee_i$estimate[[1]], 3), round(mi96.coffee_r$estimate[[1]], 3), 
           round(mi96.grass_i$estimate[[1]], 3), round(mi96.grass_r$estimate[[1]], 3), round(mi96.legum_i$estimate[[1]], 3), 
           round(mi96.legum_r$estimate[[1]], 3), round(mi96.maize_i$estimate[[1]], 3), round(mi96.maize_r$estimate[[1]], 3), 
           round(mi96.rice_i$estimate[[1]], 3), round(mi96.rice_r$estimate[[1]], 3), round(mi96.soy_i$estimate[[1]], 3), 
           round(mi96.soy_r$estimate[[1]], 3), round(mi96.sugar_i$estimate[[1]], 3), round(mi96.sugar_r$estimate[[1]], 3), 
           round(mi96.amphib$estimate[[1]], 3), round(mi96.birds$estimate[[1]], 3), round(mi96.mammal$estimate[[1]], 3),
           round(mi96.urbdist$estimate[[1]], 3), round(mi96.urbexte$estimate[[1]], 3), 
           round(mi96.br230$estimate[[1]], 3), round(mi96.mainrds$estimate[[1]], 3), round(mi96.roads$estimate[[1]], 3), 
           round(mi96.elf$estimate[[1]], 3), round(mi96.popc$estimate[[1]], 3), round(mi96.popd$estimate[[1]])),
  # P-values (1996)
  pvalue96=c(mi96.ideol$p.value, mi96.defo$p.value, mi96.rain$p.value, mi96.eva$p.value,  mi96.temp$p.value,
             mi96.climaggr$p.value, mi96.alti$p.value, mi96.slope$p.value,
             mi96.access$p.value, mi96.workab$p.value, mi96.nutri$p.value, 
             mi96.water$p.value, mi96.veget$p.value, mi96.cacao_i$p.value, 
             mi96.cacao_r$p.value, mi96.coffee_i$p.value, mi96.coffee_r$p.value, 
             mi96.grass_i$p.value, mi96.grass_r$p.value, mi96.legum_i$p.value, 
             mi96.legum_r$p.value, mi96.maize_i$p.value, mi96.maize_r$p.value, 
             mi96.rice_i$p.value, mi96.rice_r$p.value, mi96.soy_i$p.value, 
             mi96.soy_r$p.value, mi96.sugar_i$p.value, mi96.sugar_r$p.value, 
             mi96.amphib$p.value, mi96.birds$p.value, mi96.mammal$p.value,
             mi96.urbdist$p.value, mi96.urbexte$p.value, 
             mi96.br230$p.value, mi96.mainrds$p.value, mi96.roads$p.value,
             mi96.elf$p.value, mi96.popc$p.value, mi96.popd$p.value),
  # Moran's I (2000)
  test00=c(round(mi00.ideol$estimate[[1]], 3), round(mi00.defo$estimate[[1]], 3), 
           round(mi00.rain$estimate[[1]], 3), round(mi00.eva$estimate[[1]], 3), round(mi00.temp$estimate[[1]], 3),
           round(mi00.climaggr$estimate[[1]], 3), round(mi00.alti$estimate[[1]], 3), round(mi00.slope$estimate[[1]], 3), 
           round(mi00.access$estimate[[1]], 3), round(mi00.workab$estimate[[1]], 3), round(mi00.nutri$estimate[[1]], 3), 
           round(mi00.water$estimate[[1]], 3), round(mi00.veget$estimate[[1]], 3), round(mi00.cacao_i$estimate[[1]], 3), 
           round(mi00.cacao_r$estimate[[1]], 3), round(mi00.coffee_i$estimate[[1]], 3), round(mi00.coffee_r$estimate[[1]], 3), 
           round(mi00.grass_i$estimate[[1]], 3), round(mi00.grass_r$estimate[[1]], 3), round(mi00.legum_i$estimate[[1]], 3), 
           round(mi00.legum_r$estimate[[1]], 3), round(mi00.maize_i$estimate[[1]], 3), round(mi00.maize_r$estimate[[1]], 3), 
           round(mi00.rice_i$estimate[[1]], 3), round(mi00.rice_r$estimate[[1]], 3), round(mi00.soy_i$estimate[[1]], 3), 
           round(mi00.soy_r$estimate[[1]], 3), round(mi00.sugar_i$estimate[[1]], 3), round(mi00.sugar_r$estimate[[1]], 3), 
           round(mi00.amphib$estimate[[1]], 3), round(mi00.birds$estimate[[1]], 3), round(mi00.mammal$estimate[[1]], 3),
           round(mi00.urbdist$estimate[[1]], 3), round(mi00.urbexte$estimate[[1]], 3), 
           round(mi00.br230$estimate[[1]], 3), round(mi00.mainrds$estimate[[1]], 3), round(mi00.roads$estimate[[1]], 3), 
           round(mi00.elf$estimate[[1]], 3), round(mi00.popc$estimate[[1]], 3), round(mi00.popd$estimate[[1]])),
  pvalue00=c(mi00.ideol$p.value, mi00.defo$p.value, mi00.rain$p.value, mi00.eva$p.value,  mi00.temp$p.value,
             mi00.climaggr$p.value, mi00.alti$p.value, mi00.slope$p.value,
             mi00.access$p.value, mi00.workab$p.value, mi00.nutri$p.value, 
             mi00.water$p.value, mi00.veget$p.value, mi00.cacao_i$p.value, 
             mi00.cacao_r$p.value, mi00.coffee_i$p.value, mi00.coffee_r$p.value, 
             mi00.grass_i$p.value, mi00.grass_r$p.value, mi00.legum_i$p.value, 
             mi00.legum_r$p.value, mi00.maize_i$p.value, mi00.maize_r$p.value, 
             mi00.rice_i$p.value, mi00.rice_r$p.value, mi00.soy_i$p.value, 
             mi00.soy_r$p.value, mi00.sugar_i$p.value, mi00.sugar_r$p.value, 
             mi00.amphib$p.value, mi00.birds$p.value, mi00.mammal$p.value,
             mi00.urbdist$p.value, mi00.urbexte$p.value,
             mi00.br230$p.value, mi00.mainrds$p.value, mi00.roads$p.value,
             mi00.elf$p.value, mi00.popc$p.value, mi00.popd$p.value),
  test04=c(round(mi04.ideol$estimate[[1]], 3), round(mi04.defo$estimate[[1]], 3), 
           round(mi04.rain$estimate[[1]], 3), round(mi04.eva$estimate[[1]], 3), round(mi04.temp$estimate[[1]], 3),
           round(mi04.climaggr$estimate[[1]], 3), round(mi04.alti$estimate[[1]], 3), round(mi04.slope$estimate[[1]], 3), 
           round(mi04.access$estimate[[1]], 3), round(mi04.workab$estimate[[1]], 3), round(mi04.nutri$estimate[[1]], 3), 
           round(mi04.water$estimate[[1]], 3), round(mi04.veget$estimate[[1]], 3), round(mi04.cacao_i$estimate[[1]], 3), 
           round(mi04.cacao_r$estimate[[1]], 3), round(mi04.coffee_i$estimate[[1]], 3), round(mi04.coffee_r$estimate[[1]], 3), 
           round(mi04.grass_i$estimate[[1]], 3), round(mi04.grass_r$estimate[[1]], 3), round(mi04.legum_i$estimate[[1]], 3), 
           round(mi04.legum_r$estimate[[1]], 3), round(mi04.maize_i$estimate[[1]], 3), round(mi04.maize_r$estimate[[1]], 3), 
           round(mi04.rice_i$estimate[[1]], 3), round(mi04.rice_r$estimate[[1]], 3), round(mi04.soy_i$estimate[[1]], 3), 
           round(mi04.soy_r$estimate[[1]], 3), round(mi04.sugar_i$estimate[[1]], 3), round(mi04.sugar_r$estimate[[1]], 3), 
           round(mi04.amphib$estimate[[1]], 3), round(mi04.birds$estimate[[1]], 3), round(mi04.mammal$estimate[[1]], 3),
           round(mi04.urbdist$estimate[[1]], 3), round(mi04.urbexte$estimate[[1]], 3), 
           round(mi04.br230$estimate[[1]], 3), round(mi04.mainrds$estimate[[1]], 3), round(mi04.roads$estimate[[1]], 3), 
           round(mi04.elf$estimate[[1]], 3), round(mi04.popc$estimate[[1]], 3), round(mi04.popd$estimate[[1]])),
  pvalue04=c(mi04.ideol$p.value, mi04.defo$p.value, mi04.rain$p.value, mi04.eva$p.value,  mi04.temp$p.value,
             mi04.climaggr$p.value, mi04.alti$p.value, mi04.slope$p.value,
             mi04.access$p.value, mi04.workab$p.value, mi04.nutri$p.value, 
             mi04.water$p.value, mi04.veget$p.value, mi04.cacao_i$p.value, 
             mi04.cacao_r$p.value, mi04.coffee_i$p.value, mi04.coffee_r$p.value, 
             mi04.grass_i$p.value, mi04.grass_r$p.value, mi04.legum_i$p.value, 
             mi04.legum_r$p.value, mi04.maize_i$p.value, mi04.maize_r$p.value, 
             mi04.rice_i$p.value, mi04.rice_r$p.value, mi04.soy_i$p.value, 
             mi04.soy_r$p.value, mi04.sugar_i$p.value, mi04.sugar_r$p.value, 
             mi04.amphib$p.value, mi04.birds$p.value, mi04.mammal$p.value,
             mi04.urbdist$p.value, mi04.urbexte$p.value,
             mi04.br230$p.value, mi04.mainrds$p.value, mi04.roads$p.value,
             mi04.elf$p.value, mi04.popc$p.value, mi04.popd$p.value),
  test08=c(round(mi08.ideol$estimate[[1]], 3), round(mi08.defo$estimate[[1]], 3), 
           round(mi08.rain$estimate[[1]], 3), round(mi08.eva$estimate[[1]], 3), round(mi08.temp$estimate[[1]], 3),
           round(mi08.climaggr$estimate[[1]], 3), round(mi08.alti$estimate[[1]], 3), round(mi08.slope$estimate[[1]], 3), 
           round(mi08.access$estimate[[1]], 3), round(mi08.workab$estimate[[1]], 3), round(mi08.nutri$estimate[[1]], 3), 
           round(mi08.water$estimate[[1]], 3), round(mi08.veget$estimate[[1]], 3), round(mi08.cacao_i$estimate[[1]], 3), 
           round(mi08.cacao_r$estimate[[1]], 3), round(mi08.coffee_i$estimate[[1]], 3), round(mi08.coffee_r$estimate[[1]], 3), 
           round(mi08.grass_i$estimate[[1]], 3), round(mi08.grass_r$estimate[[1]], 3), round(mi08.legum_i$estimate[[1]], 3), 
           round(mi08.legum_r$estimate[[1]], 3), round(mi08.maize_i$estimate[[1]], 3), round(mi08.maize_r$estimate[[1]], 3), 
           round(mi08.rice_i$estimate[[1]], 3), round(mi08.rice_r$estimate[[1]], 3), round(mi08.soy_i$estimate[[1]], 3), 
           round(mi08.soy_r$estimate[[1]], 3), round(mi08.sugar_i$estimate[[1]], 3), round(mi08.sugar_r$estimate[[1]], 3), 
           round(mi08.amphib$estimate[[1]], 3), round(mi08.birds$estimate[[1]], 3), round(mi08.mammal$estimate[[1]], 3),
           round(mi08.urbdist$estimate[[1]], 3), round(mi08.urbexte$estimate[[1]], 3), 
           round(mi08.br230$estimate[[1]], 3), round(mi08.mainrds$estimate[[1]], 3), round(mi08.roads$estimate[[1]], 3), 
           round(mi08.elf$estimate[[1]], 3), round(mi08.popc$estimate[[1]], 3), round(mi08.popd$estimate[[1]])),
  pvalue08=c(mi08.ideol$p.value, mi08.defo$p.value, mi08.rain$p.value, mi08.eva$p.value,  mi08.temp$p.value,
             mi08.climaggr$p.value, mi08.alti$p.value, mi08.slope$p.value,
             mi08.access$p.value, mi08.workab$p.value, mi08.nutri$p.value, 
             mi08.water$p.value, mi08.veget$p.value, mi08.cacao_i$p.value, 
             mi08.cacao_r$p.value, mi08.coffee_i$p.value, mi08.coffee_r$p.value, 
             mi08.grass_i$p.value, mi08.grass_r$p.value, mi08.legum_i$p.value, 
             mi08.legum_r$p.value, mi08.maize_i$p.value, mi08.maize_r$p.value, 
             mi08.rice_i$p.value, mi08.rice_r$p.value, mi08.soy_i$p.value, 
             mi08.soy_r$p.value, mi08.sugar_i$p.value, mi08.sugar_r$p.value, 
             mi08.amphib$p.value, mi08.birds$p.value, mi08.mammal$p.value,
             mi08.urbdist$p.value, mi08.urbexte$p.value,
             mi08.br230$p.value, mi08.mainrds$p.value, mi08.roads$p.value,
             mi08.elf$p.value, mi08.popc$p.value, mi08.popd$p.value),
  row.names=c("Municipal ideological score", "Deforested area", "Rainfall", "Evapotranspiration", "Temperature",
              "Climate aggressiveness", "Altitude", "Slope",
              "Accessibility", "Workability", "Nutrients",
              "Water bodies", "Vegetation", "Cacao suit. (irrigated)",
              "Cacao suit. (rain-fed)", "Coffee suit. (irrigated)", "Coffee suit. (rain-fed)", 
              "Pasture grasses suit. (irrigated)", "Pasture grasses suit. (rain-fed)", "Pasture legumes suit. (irrigated)",
              "Pasture legumes suit. (rain-fed)", "Maize suit. (irrigated)", "Maize suit. (rain-fed)",
              "Rice suit. (irrigated)",  "Rice suit. (rain-fed)", "Soybeans suit. (irrigated)",
              "Soybeans suit. (rain-fed)", "Sugarcane suit. (irrigated)", "Sugarcane suit. (rain-fed)",
              "Threatened amphibians", "Threatened birds", "Threatened mammals",
              "Min. dist. from urban areas", "Urban", 
              "Min. dist. from Transamazonica", "Min. dist. from federal highway", "Min. dist. from federal roads",
              "Ethnolinguistic fractionalization index", "Population count", "Population density")
)


addtorow <- list()
addtorow$pos <- list(0, 0, 40)
addtorow$command <- c(paste0(paste0("& \\multicolumn{2}{c}{1996} & \\multicolumn{2}{c}{2000} & \\multicolumn{2}{c}{2004} & \\multicolumn{2}{c}{2008}", collapse=''), '\\\\'),
                      paste0(paste0("& diff. & p-val. & diff. & p-val. & diff. & p-val. & diff. & p-val.", collapse=''), '\\\\'),
                      paste0(paste0(" \\emph{N} & \\multicolumn{2}{c}{45420} & \\multicolumn{2}{c}{52918} & \\multicolumn{2}{c}{64194} & \\multicolumn{2}{c}{64664}", collapse=''), '\\\\')
)

sink("tables/mi_stats.tex")
print(xtable(mi.tab, align=c("l", rep("c",8)), digits=3),
      include.colnames = FALSE, 
      add.to.row = addtorow, 
      floating = FALSE, 
      booktabs = TRUE,
      size = "small",
      table.placement="h!"
)
sink()



### (F) FIGURES: MAIN PLOTS ----------------------------------------------------
### "Findings" section of the paper
### NOTE: first run the relevant models in 'analysis.R' 

## Figure 3: Marginal effect of district-level vote share for the incumbent’s president party on the impact of Coalition Alignment ----
## Model 3, Table A6
c.plot <- marg.effect(core[3], df.s$pres_vote_pre, c(1, 3), "p_coal") +
  xlab("Prior Vote Share")
pdf(file = "figures/marg_eff_voteshare.pdf")
ggMarginal(c.plot, type = 'histogram', margins = 'x')
dev.off()

## Figure 4: Average annual extraction by municipalities with and without protected areas before and after the commodities boom ----
# Greater than the median
df.ex$pa.median <- as.character(ifelse(df.ex$protected_prop > 0, "Yes", "No"))
df.exgr <- group_by(df.ex, year, pa.median) %>% 
  summarise(soy = mean(log(soybeantn + 1), na.rm = TRUE),
            minerals = mean(log(minleases_den + 1), na.rm = TRUE)) %>%
  filter(year != 1996)
# Parallel trends: soybean production
pt.soy <- ggplot(data = df.exgr, aes(x = year, y = soy, group = pa.median))
pt.soy <- pt.soy + geom_smooth(aes(color = pa.median)) 
pt.soy <- pt.soy + scale_x_continuous(breaks = c(1997:2010), labels = c(1997:2010))
pt.soy <- pt.soy + geom_vline(xintercept = 2001, linetype = "dashed", color = "black", size = .75)
pt.soy <- pt.soy + labs(title = " ", x = "Year", y = "Mean of Log Soybean Production", color = "Protected Areas")
pdf(file = "figures/pt_soybeans.pdf")
pt.soy
dev.off()
# Parallel trends: mining leases
pt.min <- ggplot(data = df.exgr, aes(x = year, y = minerals, group = pa.median))
pt.min <- pt.min + geom_smooth(aes(color = pa.median)) 
pt.min <- pt.min + scale_x_continuous(breaks = c(1997:2010), labels = c(1997:2010))
pt.min <- pt.min + geom_vline(xintercept = 2001, linetype = "dashed", color = "black", size = .75)
pt.min <- pt.min + labs(title = " ", x = "Year", y = "Mean of Log Mining Leases", color = "Protected Areas")
pdf(file = "figures/pt_minerals.pdf")
pt.min
dev.off()


## Figure 5: Treatment Effect for Federal, Indigenous, and State Protected Areas ----
# Save estimates ('type')
m.coef <- sapply(type, FUN = function(x) x$coefficients[1])
m.se <- sapply(type, FUN = function(x) x$cse[1])
# Data frame for model plot 
model.plot <- data.frame(model = paste(rep(c("Federal", "Indigenous", "State"), each = 3), 
                                       rep(c("Model 1", "Model 2", "Model 3"), time = 3), sep=" "),
                         coef = m.coef,
                         upper = m.coef + 1.96*m.se,
                         lower = m.coef - 1.96*m.se,
                         color = rep(c("Federal", "Indigenous", "State"), each = 3))
model.plot$model <- factor(model.plot$model, levels = model.plot$model)
# Plot
types.plot <- ggplot(data = model.plot, aes(model, coef)) +
  geom_ribbon(aes(ymin = lower, ymax = upper, fill = factor(color), color = factor(color))) +
  geom_point() + geom_point(aes(model, upper), show.legend = FALSE, size = .75, shape = "square") +
  geom_point(aes(model, lower), show.legend = FALSE, size = .75, shape = "square") +
  geom_hline(yintercept = 0)  +
  theme(axis.text.x = element_text(angle=50, vjust=0.6)) +
  theme(legend.title=element_blank()) +
  labs(x = "Model", y = "Coefficient")
pdf(file = "figures/types.pdf")
types.plot
dev.off()



### (G) FIGURES: EXPLORATORY ANALYSES ------------------------------------------

## Figure A6: Nonlinear marginal effects of Coalition Alignment at different levels of presidential vote share ----
# Save estimates (full model 3 of 'cbin')
b1.coef <- cbin[[1]]$coefficients[1]
b3.coef <- cbin[[1]]$coefficients[10:17]
b1.var <- diag(cbin[[1]]$clustervcv)[1]
b3.var <- diag(cbin[[1]]$clustervcv)[10:17]
b1b3.covar <- cbin[[1]]$clustervcv[1, 10:17]
marg <- b1.coef + b3.coef
m_se <- sqrt(b1.var + (b3.var) + (2 * b1b3.covar))
# Data frame for model plot 
model.plot <- data.frame(model = c("20-30%", "30-40%", "40-50%", "50-60%", "60-70%", "70-80%", "80-90%", "> 90%"),
                         marg = marg,
                         upper = marg + 1.96*m_se,
                         lower = marg - 1.96*m_se,
                         color = "black")
model.plot$model <- factor(model.plot$model, levels = model.plot$model)
# Plot
core.plot <- ggplot(data = model.plot, aes(model, marg)) +
  geom_segment(aes(y = lower, yend = upper, xend = model), color = "grey75") +
  geom_point() + geom_point(aes(model, upper), show.legend = FALSE, size = .75, shape = "square") +
  geom_point(aes(model, lower), show.legend = FALSE, size = .75, shape = "square") +
  geom_hline(yintercept = 0)  +
  theme(axis.text.x = element_text(angle=50, vjust=0.6)) +
  theme(legend.title=element_blank()) +
  labs(x = "Incumbent President Vote Share", y = "Estimated Marginal Effect of Alignment")
pdf(file = "figures/core_effects.pdf")
core.plot
dev.off()


## Figure A7: Marginal effect of district-level margin of victory for the incumbent’s president party on the impact of Coalition Alignment ----
## Model 3, Table A7
m.plot <- marg.effect(vmarg[3], df.s$pres_marg_pre, c(1, 3), "p_coal") +
  xlab("Prior Margin of Victory")
pdf(file = "figures/marg_eff_votemarg.pdf")
ggMarginal(m.plot, type = 'histogram', margins = 'x')
dev.off()


## Figure A8: Marginal effect of prior deforestation on the impact of Coalition Alignment ----
## Model 1, Table A9
## Deforestation interaction effect
p <- marg.effect(int[1], df.s$deforested, c(1, 4), "p_coal") + 
  xlab("Deforestation")
pdf(file = "figures/marg_eff.pdf")
ggMarginal(p, type = 'histogram', margins = 'x')
dev.off()


## Figure A9: Marginal effect of distance from Transamazônica on the impact of Coalition Alignment ----
## Model 1, Table A10
p.rds <- marg.effect(int.rds[1], df.s$br230_distance_1993, c(1, 4), "p_coal") + 
  xlab("Distance from Transamazonica")
pdf(file = "figures/marg_eff_rds.pdf")
ggMarginal(p.rds, type = 'histogram', margins = 'x')
dev.off()


## Figure A10: Marginal effect of the minimum distance from any of the Amazon’s five main federal roadways on the impact of Coalition Alignment ----
## Model 2, Table A10
p.mfr <- marg.effect(int.rds[2], df.s$main_rds_dist, c(1, 4), "p_coal") + 
  xlab("Distance from a Main Federal Highway")
pdf(file = "figures/marg_eff_mfr.pdf")
ggMarginal(p.mfr, type = 'histogram', margins = 'x')
dev.off()



### (H) FIGURES: EXTENSIONS ----------------------------------------------------

## Figure A11: Average annual peasant agriculture by municipalities with and without protected areas before and after the commodities boom ----
# Greater than the median
df.ex$pa.median <- as.character(ifelse(df.ex$protected_prop > 0, "Yes", "No"))
df.exgr <- group_by(df.ex, year, pa.median) %>% 
  summarise(families = mean(log(familiesn + 1), na.rm = TRUE),
            farmsteads = mean(log(peasantf + 1), na.rm = TRUE)) %>%
  filter(year != 1996)
# Parallel trends: peasant families
pt.fam <- ggplot(data = df.exgr, aes(x = year, y = families, group = pa.median))
pt.fam <- pt.fam + geom_smooth(aes(color = pa.median)) 
pt.fam <- pt.fam + scale_x_continuous(breaks = c(1997:2010), labels = c(1997:2010))
pt.fam <- pt.fam + geom_vline(xintercept = 2001, linetype = "dashed", color = "black", size = .75)
pt.fam <- pt.fam + labs(title = " ", x = "Year", y = "Mean of Log Peasant Families", color = "Protected Areas")
pdf(file = "figures/pt_families.pdf")
pt.fam
dev.off()
# Parallel trends: peasant farmsteads
pt.far <- ggplot(data = df.exgr, aes(x = year, y = farmsteads, group = pa.median))
pt.far <- pt.far + geom_smooth(aes(color = pa.median)) 
pt.far <- pt.far + scale_x_continuous(breaks = c(1997:2010), labels = c(1997:2010))
pt.far <- pt.far + geom_vline(xintercept = 2001, linetype = "dashed", color = "black", size = .75)
pt.far <- pt.far + labs(title = " ", x = "Year", y = "Mean of Log Peasant Farmsteads", color = "Protected Areas")
pdf(file = "figures/pt_farmsteads.pdf")
pt.far
dev.off()


## Figure A12: Marginal effect of prior municipal deforestation on the impact of Coalition Alignment ----
## Model 3, Table A25
p.md <- marg.effect(md[1], df.s$deforested_muni_pr, c(1, 4), "p_coal") + 
  xlab("Municipal Deforestation")
pdf(file = "figures/marg_eff_md.pdf")
ggMarginal(p.md, type = 'histogram', margins = 'x')
dev.off()


## Figure A13: Marginal effect of critical areas of deforestation on the impact of Coalition Alignment ----
## Model 6, Table A25
# Save estimates
b1.coef <- md[[6]]$coefficients[1]
b3.coef <- md[[6]]$coefficients[3]
b1.var <- diag(md[[6]]$clustervcv)[1]
b3.var <- diag(md[[6]]$clustervcv)[3]
b1b3.covar <- md[[6]]$clustervcv[1, 3]
marg <- b1.coef + b3.coef
b1_se <- sqrt(b1.var)
m_se <- sqrt(b1.var + (b3.var) + (2 * b1b3.covar))
# Data frame for plot
model.plot <- data.frame(model = c("No Critical Municipality", "Critical Municipality"),
                         marg = c(b1.coef, marg),
                         upper = c(c(b1.coef + 1.96*b1_se), c(marg + 1.96*m_se)),
                         lower = c(c(b1.coef - 1.96*b1_se), c(marg - 1.96*m_se)),
                         color = "black")
model.plot$model <- factor(model.plot$model, levels = model.plot$model)
# Plot
ca.plot <- ggplot(data = model.plot, aes(model, marg)) +
  geom_segment(aes(y = lower, yend = upper, xend = model), color = "grey75") +
  geom_point() + geom_point(aes(model, upper), show.legend = FALSE, size = .75, shape = "square") +
  geom_point(aes(model, lower), show.legend = FALSE, size = .75, shape = "square") +
  geom_hline(yintercept = 0)  +
  theme(axis.text.x = element_text(angle=50, vjust=0.6)) +
  theme(legend.title=element_blank()) +
  labs(x = "Critical Areas of Deforestation", y = "Estimated Marginal Effect of Alignment")
pdf(file = "figures/marg_eff_ca.pdf")
ca.plot
dev.off()



### (I) OTHER FIGURES ----------------------------------------------------------
### Identifying Assumptions
### Time Series

## Figure A3: Balance test for Coalition Alignment on 40 pre-treatment covariates in the 1996 mayoral election ----
# Assign covariate labels
vnames <- rbind("Municipal ideological score", "Deforested area", "Rainfall", "Evapotransporation", "Temperature", "Climate aggressiveness", "Altitude", "Slope",
                "Accessibility", "Workability", "Nutrients", "Water bodies", "Vegetation", 
                "Cacao suitability (irrigated)", "Cacao suitability (rain-fed)", "Coffee suitability (irrigated)",
                "Coffee suitability (rain-fed)", "Grasses suitability (irrigated)", "Grasses suitability (rain-fed)", "Legumes suitability (irrigated)",
                "Legumes suitability (rain-fed)", "Maize suitability (irrigated)", "Maize suitability (rain-fed)",
                "Rice suitability (irrigated)", "Rice suitability (rain-fed)", "Soybeans suitability (irrigated)", "Soybeans suitability (rain-fed)",
                "Sugarcane suitability (irrigated)", "Sugarcane suitability (rain-fed)", "Threatened amphibians", "Threatened birds",
                "Threatened mammals", "Minimum distance from urban areas", "Urban", "Minimum distance from Transamazonica", "Minimum distance from federal highway", "Mininum distance from federal road", 
                "Ethnolinguistic fractionalization index", "Population count","Population density")
# Means (treatment)
means.t1 <- round(rbind(
  mean(df_1996$ideo_imp[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$deforested[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$rain_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(scale(df_1996$eva_mean)[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$temp_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(scale(df_1996$climate_aggr_lv)[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$altg_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(scale(df_1996$slp_index)[df_1996$p_coal==1], na.rm=TRUE),
  mean(scale(df_1996$acc_mean)[df_1996$p_coal==1], na.rm=TRUE),
  mean(scale(df_1996$wrk_mean)[df_1996$p_coal==1], na.rm=TRUE),
  mean(scale(df_1996$nut_mean)[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$water_area[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$vegetation_area[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$cci_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$ccr_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$cfi_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$cfr_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$pgi_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$pgr_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$pli_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$plr_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$mzi_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$mzr_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$rci_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$rcr_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$sbi_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$sbr_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$sgi_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$sgr_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$amp_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$brd_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$mml_mean[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$HubDist[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$urb_extent[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$br230_distance_1993[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$main_rds_dist[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$rds_dist[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$elf[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$pop95c_avg[df_1996$p_coal==1], na.rm=TRUE),
  mean(df_1996$pop95d_avg[df_1996$p_coal==1], na.rm=TRUE)),
  digits=2)
# Means (control)
means.t0 <- round(rbind(
  mean(df_1996$ideo_imp[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$deforested[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$rain_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(scale(df_1996$eva_mean)[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$temp_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(scale(df_1996$climate_aggr_lv)[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$altg_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(scale(df_1996$slp_index)[df_1996$p_coal==0], na.rm=TRUE),
  mean(scale(df_1996$acc_mean)[df_1996$p_coal==0], na.rm=TRUE),
  mean(scale(df_1996$wrk_mean)[df_1996$p_coal==0], na.rm=TRUE),
  mean(scale(df_1996$nut_mean)[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$water_area[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$vegetation_area[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$cci_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$ccr_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$cfi_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$cfr_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$pgi_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$pgr_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$pli_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$plr_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$mzi_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$mzr_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$rci_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$rcr_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$sbi_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$sbr_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$sgi_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$sgr_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$amp_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$brd_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$mml_mean[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$HubDist[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$urb_extent[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$br230_distance_1993[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$main_rds_dist[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$rds_dist[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$elf[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$pop95c_avg[df_1996$p_coal==0], na.rm=TRUE),
  mean(df_1996$pop95d_avg[df_1996$p_coal==0], na.rm=TRUE)), 
  digits=2)
# P-values
pvals <- round(rbind(
  c96.ideol$cpval, c96.defo$cpval, c96.rain$cpval, c96.eva$cpval, c96.temp$cpval, c96.climaggr$cpval, c96.alti$cpval, 
  c96.slope$cpval, c96.access$cpval, c96.workab$cpval, c96.nutri$cpval, c96.water$cpval, 
  c96.veget$cpval, c96.cacao_i$cpval, c96.cacao_r$cpval, c96.coffee_i$cpval, c96.coffee_r$cpval, 
  c96.grass_i$cpval, c96.grass_r$cpval, c96.legum_i$cpval, c96.legum_r$cpval, c96.maize_i$cpval, 
  c96.maize_r$cpval, c96.rice_i$cpval, c96.rice_r$cpval, c96.soy_i$cpval, c96.soy_r$cpval, 
  c96.sugar_i$cpval, c96.sugar_r$cpval, c96.amphib$cpval, c96.birds$cpval, c96.mammal$cpval, 
  c96.urbdist$cpval, c96.urbexte$cpval, c96.br230$cpval, c96.mainrds$cpval, c96.roads$cpval, 
  c96.elf$cpval, c96.popc$cpval, c96.popd$cpval),
  digits=2)
# Generate plot
results <- cbind(vnames, means.t1, means.t0, pvals)
results <- results[order(as.numeric(results[, 4])), ]
plot.pval(results=results, title="", legend=FALSE, legendx=0.15, legendy=2.2, textsize=1, parcex=0.8, at1=-0.35, at2=-0.15, at3=-0.9, xlim1=-0.85)



## Figure A4: Protected Areas by Type, 1997-2012 ----
# Load cell data on protected areas
ts <- read.csv(paste0(user, "data/analysis/grid_yr dataset/protect_ts.csv"))

# Cumulative protected areas (km2)
ts.all <- data.frame(year = rep(ts$year, each = 3),
                     type = rep(c("fed", "ind", "st"), times = length(ts$year)),
                     prot = NA,
                     stringsAsFactors = FALSE)
for (i in 1:nrow(ts.all)){
  selc <- which(colnames(ts) == ts.all$type[i])
  selr <- which(ts$year == ts.all$year[i])
  ts.all$prot[i] <- ts[selr, selc]
}
ts.all$type <- rep(c("Federal", "Indigenous", "State"), times = length(ts$year))

# New protected areas (km2)
ts.new <- data.frame(year = rep(ts$year, each = 3),
                     type = rep(c("fed.new", "ind.new", "st.new"), times = length(ts$year)),
                     prot = NA,
                     stringsAsFactors = FALSE)
for (i in 1:nrow(ts.new)){
  selc <- which(colnames(ts) == ts.new$type[i])
  selr <- which(ts$year == ts.new$year[i])
  ts.new$prot[i] <- ts[selr, selc]
}
ts.new$type <- rep(c("Fed.", "Ind.", "State"), times = length(ts$year))
ts.new$prot <- round(ts.new$prot)

# Plot: cumulative areas
pl.all <- ggplot(ts.all, aes(year, prot)) + geom_col(aes(fill = type), width = 0.4) + 
  scale_x_continuous(breaks=c(1997:2012), labels=c(1997:2012), limits=c(1996.5, 2012.5)) +
  ylab("Protected Area (Km2)") + xlab("Year") + theme(legend.position = "none") + 
  scale_y_continuous(labels = comma)
pdf(file = "figures/ts_all.pdf")
pl.all
dev.off()
# Plot: new areas
pl.new <- ggplot(ts.new, aes(year, prot)) + geom_col(aes(fill = type), width = 0.4) + 
  scale_x_continuous(breaks=c(1997:2012), labels=c(1997:2012), limits=c(1996.5, 2012.5)) +
  xlab("Year") + theme(legend.title = element_blank(), axis.title.y = element_blank()) + 
  geom_vline(xintercept = c(1998.5, 2002.5, 2006.5, 2010.5), color = "blue", linetype = 2) + 
  geom_vline(xintercept = c(2000.5, 2004.5, 2008.5), color = "red", linetype = 6) + 
  scale_y_continuous(labels = comma)
pdf(file = "figures/ts_new.pdf")
pl.new
dev.off()


## Figure A5: Time-series shifts in coalition alignment between president and mayors, 1997-2012 ----
# Load municipality data on alignments
alignm <- read.csv(paste0(user, "data/analysis/grid_yr dataset/coalitions.csv"))
# Reshape
alignm <- reshape(alignm, varying = c(2:29), direction = "long", idvar = "codigo", sep = ".", timevar = "year") 
# Coalitional changes groups
gov.opp <- group_by(alignm, year) %>% 
  summarise(total = n(), 
            govt = sum(stgov), 
            oppt = sum(stopp), 
            togovt = sum(togov), 
            tooppt = sum(toopp))
# Time series: frequencies and proportions
ts.qnt <- rbind(as.data.frame(t(gov.opp[1,3:6])), as.data.frame(t(gov.opp[2,3:6])), 
                as.data.frame(t(gov.opp[3,3:6])), as.data.frame(t(gov.opp[4,3:6])), 
                as.data.frame(t(gov.opp[5,3:6])), as.data.frame(t(gov.opp[6,3:6])),
                as.data.frame(t(gov.opp[7,3:6])))
ts.alignm <- data.frame(year = rep(gov.opp$year, each = 4),
                       type = rep(c("Stayed Gov't", "Stayed Opp", "To Gov't", "To Opp"), times = length(gov.opp$year)),
                       stringsAsFactors = FALSE) %>% cbind(ts.qnt)
colnames(ts.alignm)[3] <- "qnt"
ts.alignm$shr <- round((ts.alignm$qnt / 666), 3)
# Plot
ts.coal <- ggplot(data=ts.alignm, aes(x=year, y=shr, fill=type)) +
  geom_bar(stat="identity", width=.75) +
  geom_vline(xintercept = c(1998.5, 2002.5, 2006.5, 2010.5), color = "blue", linetype = 2) + 
  geom_vline(xintercept = c(2000.5, 2004.5, 2008.5), color = "red", linetype = 6) + 
  scale_x_continuous(breaks=c(seq(1998,2010,1)), labels=c(seq(1998,2010,1))) +
  scale_fill_manual(values=c("#999999", "#E69F00", "#56B4E9", "#009E73")) +
  labs(x = "Year", y = "Share of municipal governments", fill = "Coalitional changes")  
pdf(file = "figures/ts_pcoal_stack.pdf")
ts.coal
dev.off()


